Configuration space connectivity across the fragile to strong transition in silica 
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We present a numerical analysis for SiC>2 of the fraction of 
diffusive direction /diff f° r temperatures T on both sides of the 
fragile-to-strong crossover. The T-dependence of /^iff clearly 
reveals this change in dynamical behavior. We find that for T 
above the crossover (fragile region) the system is always close 
to ridges of the potential energy surface (PES), while below 
the crossover (strong region), the system mostly explores the 
PES local minima. Despite this difference, the power law 
dependence of f^g on the diffusion constant, as well as the 
power law dependence of f^g on the configurational entropy, 
shows no change at the fragile to strong crossover. 



PACS numbers: 61.43.Fs, 64.70.Pf, 66.10.Cb 



A useful approach for relating dynamics and thermo- 
dynamics of glassforming liquids in their supercooled 
states is offered by the study of the liquid's potential en- 
ergy surface (PES) 0. The PES is the multidimensional 
surface generated by the system's potential energy as a 
function of all atomic coordinates. Numerical and theo- 
retical studies P-[l0| are providing evidence that several 
aspects of the liquid dynamics have a clear signature in 
the properties of the explored regions of the PES. For ex- 
ample, correlation functions display stretching in time in 
the same temperature (T) range in which the system ex- 
plores regions of the PES associated with local minima of 
deeper and deeper energy |flj . Similarly, the rapid slow- 
ing down of the dynamics which takes place in the weakly 
supercooled region is associated to the exploration of re- 
gions of the PES close to saddles of lower and lower order 
[P~2|-]l7|| . Recent numerical studies for fragile liquids |is|j 
have suggested the possibility that a cross over T, T x 
marks a change in the geometrical properties of the PES 
explored. In this picture, above T x the system trajecto- 
ries are mostly located close to the ridges between differ- 
ent PES basins ("border dynamics"), while below T x the 
system only rarely samples the regions connecting differ- 
ent PES basins ( "minimum-to- minimum dynamics" ) . T x 
has been associated |ll||(],|l3],|l4| with the critical tem- 
perature predicted by the ideal mode-coupling theory 
(MCT) g|. 

Instantaneous normal mode (INM) analysis jl2j is a 
powerful technique to investigate the connectivity prop- 
erties of the PES explored at different temperatures. The 



curvature of the PES along each of the 3 AT — 3 indepen- 
dent directions (the eigenvectors of the Hessian matrix) 
is calculated and analyzed to estimate the fraction of dif- 
fusive directions /dig- If the representative configuration 
is crossing a ridge separating different PES basins, some 
of the local curvatures are negative. A computationally- 
expensive but well-defined screening procedure to sort 
out the negative-curvature directions contributing to dif- 
fusion has been developed 22 2^,|ljl. For all model 
potentials for which such analysis has been performed, 
it has been shown that the T and density p variation 
of /diff controls the behavior of the long time dynamics, 
supporting the hypothesis that information about the lo- 
cal properties of the PES may be sufficient to describe 
long-time dynamical processes. 

The analysis of the PES connectivity has been limited 
to models for fragile glass-forming liquids and to tem- 
peratures above T x , due to the difficulty in generating 
equilibrium configurations at low T. In all studied cases, 
it was consistently found that T x locates the T at which 
/diff appears to extrapolate to zero Jl3| , p3|p^ ]. 

Here, we present an evaluation of the fraction of dif- 
fusive directions (using the INM analysis) for the well- 
studied BKS model for silica j24]], for which we generate 
equilibrium configurations for temperature both above 
and below T x . The INM spectrum of BKS silica has 
been previously calculated by Benbenek and Laird |2J| 
for equilibrium states above T x but no evaluation of /diff 
was reported. For the BKS model, the high T dynam- 
ics has been shown to be consistent with the predictions 
of MCT at both a qualitative and quantitative p7f| 
level. At lower T, the T-dependence of the characteris- 
tic times |p6| shows a crossover toward an Arrhenius T 
dependence which has been interpreted as a clear case of 
fragile-to-strong transition |28|j . The crossover temper- 
ature T x , around 3330 K at p = 2.36 g/cm 3 , has been 
interpreted by Horbach et ai at as the MCT critical tem- 
perature |26|. In the Arrhenius region, the activation 
energy is 54000AT and 60000AT , for oxygen and silicon 
respectively. 

We find that both the T dependence of the diffusion 
constant D and of /diff show a clear signature of the two 
different dynamical behaviors — fragile above and strong 
below T x — . Despite these differences, the relation be- 
tween D and /diff is not sensitive to the presence of dy- 
namical changes. Moreover, we find that /diff does not 
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vanish around T x but changes to an Arrhenius T depen- 
dence, similar to the T dependence of D. Finally, we 
show that the number of diffusive directions is related to 
the number of basins of the PES, providing a possible 
explanation for the recently-observed validity j|,^8|,^9| of 
the Adam-Gibbs equation ]3(| in the T region where bor- 
der dynamics dominates. 

Our results are based on extensive simulations of a 
system of 999 atoms, for p — 2.36 g/cm 3 , close to the 
density of ambient pressure silica. We investigate eleven 
T, in the range from T = 2650Jf to T = 7QQQK. We 
carry out simulations in the NVE ensemble, using a 1 fs 
integration timestep. We evaluate the long range interac- 
tion by implementing the Ewald summation. To guaran- 
tee proper equilibrium conditions all simulations lasted 
longer than several times the slowest collective structural 
relaxation time; low T runs lasted longer than 50 ns. We 
also performed averages over eight different realizations. 
For each studied state point we calculate eigenvalues and 
eigenvectors for 96 configurations. For each of the eleven 
T studied, we performed approximately 2600 minimiza- 
tions, using a total of about 50,000 CPU hours. 

The INM analysis requires the evaluation of the eigen- 
vectors and associated eigenvalues of the potential en- 
ergy second derivative matrix, the Hessian. According 
to the procedure described in Refs. |2^,||l|], all eigenvec- 
tors associated with negative eigenvalues are inspected 
in order to eliminate those associated with intra-basin 
anharmonicities ("shoulder" and "false barrier" modes). 
Fig. [j] shows (i) the T dependence of the fraction of di- 
rections with negative eigenmodes / u ,(ii) the fraction of 
directions whose one-dimensional profile is double-well 
shaped /d w , and (iii) the number of diffusive directions 
/diff, calculated by elimina ting from the double-well set 
all false barrier modes HHf. Above T = 3330K, all 
three quantities show a fast decrease with T. While 
f u and /dw assume non-zero values, /diff appears to ap- 
proach zero on cooling. A clear change of concavity in the 
T dependence of /diff takes place above T = 3330i"T. As 
found previously for fragile liquids for T > T x , the fast 
decrease of /diff confirms that above T = 333QK: (i) the 
slowing down of dynamics is associated with a progressive 
decrease of the number of possible directions that lead to 
a different basin; (ii) the system trajectories are located 
close to PES ridges (border dynamics); (iii) even for BKS 
silica, the dynamics properly described by the ideal MCT 
are border dynamics, supporting the identification of T x 
with the MCT critical temperature. 

For T below 3330K, the system spends most of the 
time far from the PES ridges. To support this state- 
ment, Fig. ^ shows the probability P(fdis) of finding a 
configurations with a specific /diff value for T = 2800 K 
(< T x ) and T = 4000 K (> T x ). Above T x all examined 
configurations have a nonzero /diff, while below T x the 
distribution is peaked around zero and most of the con- 
figurations are characterized by the absence of diffusive 



directions. 

We next investigate the functional relation between 
/diff and D. The T dependence of both quantities, shown 
in an Arrhenius plot in Figs. 0(a) and ||(b) , show a change 
in the T-dependence above and below T x . Despite this 
different dynamical behavior, the quantity /diff is related 
to D by the same power law relation both for T > T x and 
for T < T x (Fig.|(c)). In the entire studied T region, D 
follows the law 

^/T~(/diffr, (i) 

with a 1.3 ± 0.2 over more than two decades in /diff 
and more then three decades in D/T. The same func- 
tional form describes the relationship between D and 
/diff both above and below T x showing that while the 
T dependence of both D and /diff is sensitive to the 
microscopic mechanisms controlling the dynamics, the 
fragile-to-strong transition does not affect in the rela- 
tion between D and /diff- We stress also that the same 
functional form has been found to describe the relation 
between D and /diff in the case of the extended simple 
point charge (SPC/E) model for water |Q, for which 
calculations were limited to T > T x . 

The interpretation of the dynamics in terms of the frac- 
tion of diffusive directions is complementary to the anal- 
ysis which attributes the slowing down of the dynamics 
to the decrease of the liquid configurational entropy S con f 
pot - In the case of BKS silica, S CO nf has been recently 
calculated and shown to describe, via the Adam-Gibbs 
relation, the slowing down of the dynamics both above 
and below T x |2?| . Since S'conf is a measure of the number 
of distinct PES basins explored by the system — if the 
description in terms of S'conf and the description in terms 
°f /diff are both valid — a relation must exist between the 
number of basins and the number of directions connecting 
them. Fig.^ shows the relation between S con f and /diff for 
the case of BKS silica. Within the numerical uncertainty, 
the fraction of diffusive directions appears to be propor- 
tional to the number of explored basins, Q = e Sconf / k s ; 
irrespective of the strong or fragile character of the dy- 
namics. Data in Fig.||are consistent with similar findings 
which were limited to the SPC/E model for water to the 
T-region above T x . The present results suggest that the 
linear relation between log /diff and S con f is not model 
dependent; indeed it has been recently derived within 
the random energy model |3jJ. The fact that this rela- 
tion holds even above T x is particularly interesting; since 
above T x the system dynamics is a dynamics of borders, 
there is no clear reason why such border dynamics should 
be well described by the Adam-Gibbs relation which fo- 
cuses on the number of basins explored as a function of 
T. The observed relation between log /diff and S con f may 
offer a key for the resolution of this apparent paradox. It 
is a challenge for future studies to find out if /diff is a po- 
tentially richer quantity for describing dynamics in deep 
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supercooled states, as the results reported here seem to 
suggest. 
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To separate the double-well directions into diffusive 

modes (modes connecting distinct minima, or inher- 
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ent structures), and false modes (intra-basin modes) we 
perform a conjugate gradient minimization with toler- 
ance 10~ 11 kJ/mol, starting from the two one dimen- 
sional double-well minima. We then calculate the dis- 
tance between the two inherent structures d = [J^C 7 *? ~~ 
r\) 2 /N] 1 / 2 , where rf is the i th atom position in inher- 
ent structure a . The distribution of d found is bimodal, 
and the two contributions are separated by more than 
one order of magnitude. Modes which contribute to the 
peak located at "large" d ~ 0.01 nm are classified as 
diffusive modes (/diff)- All other modes are classified as 
"false" double-well modes, since both sides of the one 
dimensional profile end up in the same minimum. The 
distribution of distances is not centered at zero, but cen- 
tered around 10~ 5 nm, due to numerical uncertainties 
in the exact location of the local minimum during the 
minimization procedure. For further information on the 
procedure, and for a critical discussion of the results ob- 
tained with this procedure, see Ref. [ pl[ . 
T. Keyes, Phys. Rev. E 62, 7905 (2000). 
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FIG. 1. Temperature dependence for the BKS model of the 
fraction of modes with negative curvature /„, the fraction of 
modes with double- well shaped one-dimensional profile /d W! 
and the fraction of diffusive directions /diff- The arrow marks 
T x for the isochore studied [1261 . 
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FIG. 2. Probability P(fdis) of finding a configuration 
with a given /dig value above T x (T = 4000K) and below T x 
(T = 2800 Ji"). While for T > T x all examined configurations 
have a nonzero /diff, below T x , the distribution is peaked 
around the origin. 
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FIG. 3. Arrhenius plot of (a) the diffusion constant D for 
Si atoms; (b) /diff- The line has slope 54, OOOA' in (a) and 
36, OOO-ft' in (b) . Part (c) shows the parametric relation D/T 
vs /dig in a log-log scale. The data are perfectly smooth 
through the transition at T x . 
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FIG. 4. Parametric plot of the fraction of diffusive modes 
/diff as a function of the configurational entropy 5 con f for the 
BKS model of silica 
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